EPJ manuscript No. 

(will be inserted by the editor) 



A Lattice Gas Coupled to Two Thermal Reservoirs: Monte Carlo 
and Field Theoretic Studies 

E. L. Praestgaard 1 , B. Schmittmann 2 , and R. K. P. Zia 2 ' 3 

1 Institute for Life Sciences and Chemistry, Roskilde University, 4000 Roskilde, Denmark 

2 Center for Stochastic Processes in Science and Engineering, Department of Physics, 
Virginia Polytechnic Institute and State University, Blacksburg, VA 24061-0435, USA 

3 Fachbereich Physik, Universitat - Gesamthochschule Essen, D-45117 Essen, Germany. 



February 1, 2008 



Abstract. We investigate the collective behavior of an Ising lattice gas, driven to non-equilibrium steady 
states by being coupled to two thermal baths. Monte Carlo methods are applied to a two-dimensional 
system in which one of the baths is fixed at infinite temperature. Both generic long range correlations 
in the disordered state and critical poperties near the second order transition are measured. Anisotropic 
scaling, a key feature near criticality, is used to extract T c and some critical exponents. On the theoretical 
front, a continuum theory, in the spirit of Landau-Ginzburg, is presented. Being a renormalizable theory, 
its predictions can be computed by standard methods of e-expansions and found to be consistent with 
simulation data. In particular, the critical behavior of this system belongs to a universality class which is 
quite different from the uniformly driven Ising model. 

PACS. 64.60.Ht Dynamic critical phenomena - 64.60.Ak Renormalization-group, fractal, and percolation 
studies of phase transitions - 05.10.Ln Monte Carlo methods - 05.70.Ln Non-equilibrium thermodynamics, 
irreversible processes. 



1 Introduction 

In recent years, there has been considerable interest in 
. the collective behavior of many body systems which, de- 
spite being far from equilibrium, have settled in time- 
independent steady states. Unlike their equilibrium coun- 
terparts, there is no simple Boltzmann factor which pro- 
vides, in general, the distribution for these non-equilibrium 
steady states, so that most of the progress to date is 
made by studying specific cases. One distinguishing fea- 
ture of these systems is their coupling to more than one 
reservoir of energy, so that there is a constant "through- 
flua?' of energy. Examples include particles driven by ex- 
ternal fields such as gravity or DC electricity, gaining en- 
ergy from these fields and losing energy thermally to their 
environment. An important class is driven diffusive sys- 
tems which were motivated physically by the prop- 
erties of fast ionic conductors B and theoretically by their 
simplicity. In particular, they are arguably the simplest 
generalization of the well-known Ising model Q to non- 
equilibrium conditions. In this paper, we investigate a 
lattice gas in which particle hopping rates are controlled 
by two thermal baths, with temperatures T' and T, first 
introduced in Q. Of course, systems subjected to two 
(or more) baths are extremely common, from food be- 
ing cooked on a stove to the earth as a whole. In most 
cases, temperature gradients (in space or time or both) are 



present, so that the "steady" states are typically inhomo- 
geneous, controlled by "external" parameters with macro- 
scopic length scales. An excellent example is Rayleigh- 
Benard cells, where macroscopic length scales are intro- 
duced through gravity as well as temperature gradients. 
Not surprisingly, large scale properties in such systems 
can be quite complex, while universal features may be well 
hidden. In contrast, this two-temperature model is macro- 
scopically homogeneous in both space and time, so that 
crucial features associated with non-equilibrium steady 
states are prominently displayed in the foreground. 

Of the many large scale (in both space and time) phe- 
nomena displayed by this class of non-equilibrium sys- 
tems, two have received considerable attention. One is 
the presence of long range correlations at all temperatures 
above criticality, despite the absence of long range inter- 
actions or dynamics (e.g., jumps) at the microscopic level. 
First revealed in the uniformly driven lattice gas |6|, they 
have also been observed in the randomly driven lattice gas 
and the two-temperature models ||] . The other remark- 
able behavior is associated with the critical point, where 
significant deviations from the Ising universality class were 
found Here, we focus on the simple 

two-temperature Ising lattice gas and present an exten- 
sive Monte Carlo study, with many results beyond those 
briefly reported in fl^| . In this model, particle-hole pairs 
in the x-direction exchange randomly, simulating coupling 
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to a bath with T' = oo. Exchanges in the y-direction do 
depend on the energy change in a manner consistent with 
coupling to a reservoir at T. As we vary T, we observe a 
second order phase transition, similar to the one found by 
Onsager Q , but located at an approximately 40% higher 
value. A careful study of the two-point correlations in the 
disordered phase verifies the long-range nature (r~ 2 ) to 
a much better extent than before H. These long-range 
correlations are "generic" ||[l6| and can be easily under- 
stood [ p"7| through a Langevin equation for a continuum 
field, which is simply a non-equilibrium version of time- 
dependent Landau-Ginzburg theory. Based on this theory, 
a renormalization group (RG) analysis can be applied to 
extract the critical behavior. Reported briefly in the past 
||l8|,[l9), the details of such an analysis are given here. Sim- 
ilar to the uniformly driven lattice gas, momenta in the 
x- and y-directions scale with different power laws. Unlike 
the uniformly driven lattice gas, where the powers differ 
by a factor of three (in two dimensions), the difference 
here is much closer to the "mean- field" value of two. Re- 
lying on these results, we carried out extensive simulations 
using a series of rectangular samples (L x ~ L 2 ), so as to 
simplify the finite size scaling analysis. We measured a 
cumulant ratio, the "magnetization" , the "specific heat" 
and "energy" fluctuations. A consistent picture emerges, 
in general agreement with theoretical predictions. In con- 
trast, using the customary square samples (L x = L y ), 
no consistent set of parameters can be found to collapse 
the data. Undoubtedly, a new scaling variable is playing 
a dominant role here. Our conclusion is that anisotropic 
scaling is crucial and that the simulations support the field 
theoretic renormalization group analysis of this system. 

The paper is organized as follows. In the next section, 
the lattice model and its continuum partner will be de- 
scribed. Simulation methods and results, as well as com- 
parisons with the theoretical predictions, may be found 
in Section 3. For completeness, we present details of the 
field theoretic renormalization group analysis (Section 4). 
A concluding section is devoted to a summary and an out- 
look. 



2 The Discrete Model and A Continuum 
Theory 

2.1 A Two Temperature Lattice Gas 



The constituents of our model are identical to the two- 
dimensional Ising lattice gas On a square lattice 
with fully periodic boundary conditions, the sites i = 
(ixjiy) may be empty or occupied by a single particle, 
so that a configuration, C, is completely specified by the 
set of occupation numbers {ni}, with each n\ being or 
1. In general, our system will be rectangular, i.e., 



with L x 7^ L y typically. The particles interact with nearest 
neighbor attraction, so that, in the Hamiltonian 



H[C] = -4J J2 



(2.1) 



<ij> 



J is positive. We choose 4 J here so that the effective part 
of H in the spin language is just — J s i s j > since s — 2n — 
1 = ±1. To simulate diffusing particles, we will allow only 
particle hops to empty nearest neighbor sites, a dynamics 
corresponding to Kawasaki spin exchange pi] . Thus, the 
total particle number is conserved and, in all cases studies 
here, fixed to be N/2, where 



N = L T L„ 



(2.2) 



Such a half-filled lattice is equivalent to a spin system 
with zero total magnetization (M — 0) , so that the critical 
point can be accessed. 

If this system is in contact with a single thermal bath 
at temperature T, then the probability for finding it in 
any configuration is given by the Boltzmann factor 



P eq [C] cxexp{-ft[C]/fc B T} 



(2.3) 



iy G [1, Ly] 



and many of its properties are well known p2| . For tem- 
peratures above the Onsager |]l5) value, T Q = 2.2692 J/fcs, 
the system is in a homogeneous, disordered state. Particle- 
particle correlations decay exponentially, governed by a 
finite correlation length £. As T — > T , £ — > oo while cor- 
relations decay with a power law. Thermodynamic quan- 
tities develop singularities as a function of various param- 
eters. Due to the constraint M — 0, the ordered state 
for T < T Q is inhomogeneous, consisting of a particle-rich 
(mainly s > 0) region co-existing with a hole-rich one 
(mainly s < 0). Given the boundary conditions, each of 
these regions is a single strip, of width L > /2, spanning the 
lattice along L < (where L > is the longer dimension, etc.). 
In this manner, the interfacial energy is minimized. Thus, 
for example, states with vertical strips will dominate in a 
system with L x > L y . Finally, the average s within one 
of these regions will be, in the thermodynamic limit, the 
spontaneous magnetization: ±m. 

Our interest in this model lies in its behavior when 
placed in contact with two baths, at different tempera- 
tures. In general, we can expect energy to flow through our 
system, from one bath to another. In the limit that the 
baths are infinitely larger than the lattice gas, this energy 
flux will settle down to a constant (on the average) and our 
system reaches a non- equ ilibrium time-independent state. 
However, unlike Eqn (2J5) above, the steady state proba- 
bility distribution, P SS [C] is not known in general. Indeed, 
P ss [C] is expected to depend on the details of the dynam- 
ics, namely, how our system is coupled to the two baths. 
Nevertheless, within a class of models, we believe that 
there are many universal properties which do not depend 
on such details. Here, we will focus on a particular model, 
deferring a discussion of the universality class until Section 
4. 

To completely specify our model, we list the rules of 
time evolution: 
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— choose a random nearest neighbor particle-hole pair 

— if the pair lies along the ^-direction, exchange them 

— if the pair lies along the y-direction, exchange them 
with probability min[l, e~ 4W / fcsT ], where ATt is the 
energy change due to the exchange. 

Note that the second rule can be replaced by one sim- 
ilar to the third, with T", the temperature of the second 
bath. Here, we have further simplified the model by setting 
T" = oo. Of course, this model may be tied "continuously" 
to the equilibrium case by lowering T' down to T. 

In our Monte Carlo simulation, a random pair is first 
chosen. If this is a particle-hole pair, the rules above apply; 
otherwise, another pair will be randomly chosen. A Monte 
Carlo step (MCS) is defined by L x L y attempts. Starting 
from, typically, a random initial configuration, the system 
is evolved for up to 5 x 10 6 MCS. Allowing the first 10 5 
MCS for the system to come to steady state, we measure 
various quantities every 10 - 100 MCS. The results quoted 
in the next section are obtained by averaging these mea- 
surements. 

As will be discussed in more detail, the second order 
phase transition is found to survive. However, the critical 
temperature is raised by approximately 40%! In addition, 
the two-particle correlation function develops power law 
decays at all temperatures in the disordered phase, while 
the critical properties are modified from the equilibrium 
Ising class 2^| . To understand such collective behavior 
in the long-time and large-scale limit, we often rely on 
continuum descriptions. This is the topic of the next sub- 
section. 



2.2 Formulation of a Langevin Equation 

In this subsection, we formulate a continuum approach, 
based on field theory. For the convenience of readers who 
may not be familiar with this approach, we first sum- 
marize the well-established steps for the Ising model in 
equilibrium. In that case, a series of continuum theories 
yield better and deeper understanding of its properties. 
The simplest continuum description starts with Landau's 
macroscopic free energy function for the magnetization M 



A(M) = -M 2 



4!' 



-M 



which focuses only on spatially homogeneous features. To 
capture spatial inhomogeneities and thermal fluctuations, 
we rely on the "mesoscopic" Landau-Ginzburg Hamilto- 
nian for a magnetization </?(x) which is a coarse-grained 
version of Si from the lattice model. Furthermore, from 
rcnormalization group studies, we learned that it is im- 
portant to formulate such theories in general dimension d 
p3f (even though our simulations are restricted to d = 2). 
Thus, we write 



H LG = J d d x {\{V v ? + T -^ 2 + 1^}- 



(2.4) 



Far from the critical point, which is modeled by vanish- 
ing r here, good agreement with data can be gleaned from 



Klg by inserting Eqn (2.4) into Eqn (2.2), i.e., P eg [i^(x)] cx 
exp {— 7Ylg}? an d relying on simple approximation schemes 
in subsequent computations. Near the phase transition, 
however, more powerful renormalization group techniques 
pi[ are needed, providing excellent predictions of long- 
wavelength properties, such as critical behavior. Further, 
these appro aches can be generalized to dynamic phenom- 
ena [E5Lpq] near equilibrium, starting with a Langevin 
equation for tp(x,t). 

For later reference, we briefly review this procedure for 
the Ising model, in the presence of a conservation law on 
the total magnetization J d d x <p(x, t). Then, the Langevin 
equation takes the form of a continuity equation: 



d t Lf(x,t) = -V- J(x,t) 



(2.5) 



The form of the current J is postulated, guided by the 
symmetries and key physical features of the microscopic 
rates. It typically consists of a deterministic term, captur- 
ing the thermodynamic forces acting on <p, and a Gaussian 
noise term which reflects the effect of thermal fluctuations. 
For the Ising model, the deterministic part ensures that <p 
relaxes towards configurations of lower energy. Moreover, 
J must be chosen such that the known equilibrium distri- 
bution, P eq [ip(x)], controls the static properties of tp(x.,t). 
This constraint forces the Langevin equation to satisfy the 
fluctuation-dissipation theorem (FDT) ]27j . Let us briefly 
recall the resulting structure. The most general form of 
the Langevin equation is 



d t p = W(<p,Vp...) + ( 



(2.6) 



Here, F is a functional of tp and its derivatives, constructed 
on physical grounds in the spirit of a Landau expansion. It 
reflects the deterministic part of the dynamics. £ denotes 
the noise term, with correlations 



(C(x, f)C(x', t')) = 2N05(x - - 



(2.7) 



is just a positive constant which measures the strength 
of the correlations. More importantly, the noise "matrix" 
N carries information about conservation laws or internal 
symmetries which are obeyed by the dynamics. For our 
purposes, it is sufficient to consider only a scalar order 
parameter tp and noise matrices which do not depend on 
<p. (More general cases are discussed in, e.g., Q and [p9|.) 
Given this basic form, a steady-state solution for the con- 
figurational probability is easily found [^8| provided F is 
"Hamiltonian" , i.e., if it can be written as N acting on the 
functional derivative of an appropriate Hamiltonian TL: 



IT 



5H 

5tp 



In this case, 
distribution, 
We will refer 
§, as "FDT 
0, plays the 
the scale for 
inition of 7i. 



the steady state is simply the equilibrium 
exp(—7i./0), irrespective of the choice of N. 
to such a dynamics, in an operational sense 
-satisfying". Note that the "noise strength," 
role of a temperature here. Since it just sets 
the energy, it may be absorbed into the def- 
In the Ising case, the microscopic dynamics 
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is isotropic (at least near the critical point) and order- 
parameter conserving, so that N oc —V 2 . Our preceding 
discussion suggests the equation of motion 

a t ^(x,t) = AV 2 ^^ + C(x,i) (2.9) 

with 



(C(x,t)) 

(C(x,t)c(x',0) 



-2AV 2 <5 (x - x') S (t - t') (2.10) 



The current is therefore J(x, t) = — AV (<5Wlg/<V)+ 7 7( x ) t), 
where the second term is the noisy part, related to ( sim- 
ply by e( x ;0 = — V ' ))(x, t). The coefficient A sets the 
time scale. In fact, symmetry considerations alone would 
have forced us into this form, e ven without invoking the 
FDT: the deterministic part of (Q, V 2 (SHlg/5(p) just 
contains the leading terms of a general Landau expansion, 
in powers of (p and V^j which are consistent with all sym- 
metries of the Ising model. 

This Langevin equation, known as Model B |?0|, de- 
scribes the critical dynamics of a conserved scalar order 
parameter. The time- independent, or static, aspects of its 
critical behavior are controlled by Hlg an d fall into the 
equilibrium Ising class. In addition, one obtains time de- 
pendent properties as well, such as the dynamic critical 
exponent z and the scaling behavior of dynamic correla- 
tion and response functions. These are characteristic for 
the Model B class. We note, for later reference, that the 
critical parameter, r, now plays the role of a diffusion 
coefficient: As T drops below T c , r becomes negative indi- 
cating that "antidiffusion" , i.e., phase segregation, occurs. 

Following these lines, we formulate a continuum the- 
ory for our non- equilibrium Ising lattice gas. Once again, 
we begin with a continuity equation for </?(x, t), so that 
the conservation law is guaranteed. Before writing an ex- 
pression for the current, however, it is paramount that 
we should summarize the symmetries of our microscopic 
theory. Like the equilibrium Ising model, it exhibits full 
translation and reflection invariance, as well as the charac- 
teristic Ising "up-down" symmetry. However, in contrast 
to its equilibrium counterpart, even near criticality we 
can no longer expect full rotational symmetry, since dis- 
tinct temperatures, T a , control exchanges along different 
directions a — 1,2, ...,d. In the simulations, d — 2 and 
Ti (= T x ) = V = oo while T 2 (= T y ) = T < oc. Clearly, 
we enjoy some freedom here in extending this model to 
general dimensions: in principle, we could imagine cou- 
plings to d different temperature baths, with temperatures 
T\ > T-2 > ... > Td- However, the most interesting and 
fundamental case remains that of just two baths, with 
T 2 =T 3 = ... = T d = T and T x = T > T. In other words, 
we select a one- dimensional subspace to be coupled to the 
higher temperature (T" = oo in simulations). The full ro- 
tational symmetry, O(d), of the Ising model is obviously 
reduced toanO(d-l) symmetry in the remaining (d— 1) 
dimensions. We will refer to the subspace associated with 
the higher temperature as "parallel" (e.g., Tj| = T") and 
the complementary subspace as "transverse" (Tj_ = T) . In 



the section on field theoretic studies, we will return briefly 
to the general case, in order to show that all nontrivial fea- 
tures are indeed captured by the simpler two-temperature 
theory. 

Given the restricted rotational invariance of our model, 
we can no longer expect the same Landau expansion for 
the currents in the parallel and the transverse subspaces. 
In particular, all V operators should be split into parallel 
(d) and transverse (Vj_) parts, accompanied by different 
coefficients. Note, however, that invariance under reflec- 
tions guarantees that gradients still appear in pairs: 
or d 2 . Thus, rV 2 "splits" into t±V 2 ± + T||<9 2 , mV 2 becomes 
m_lV 2 -|-M||9 2 , and V 4 turns into a±S7 4 L -\-2a x d 2 S7 2 L +a\\d i . 
By appropriate rescaling of both parallel and transverse 
lengths, a± and aii can be set to unity. Thus, we can write 
the most general Landau expansion which still respects all 
symmetries of the microscopic theory as 

dMx,t) = -V- J(x,t) 

= A{(r ± -V 2 L )V 2 L ^+(T||-d 2 )dV 

-2 ax a 2 Vi^+^(Vi+a 3 9 2 V 3 } 



-V-tj(x,t) 



(2.11) 



where, for later convenience, we have written {u, ua 3 ) for 
(u±, Mm). Turning to the noise, while keeping it Gaussian, 
we should expect different variances in the two subspaces. 
Hence, 

(77 Q (x,^(x',i , )) = 2Aa a ^ /3 <5(x-x')^(t-i') (2.12) 

where the cr's represent the noise strengths. Since the 
transverse subspace is still fully isotropic, we have o 2 = 
... = o d = o j_ . Generically, however, we should expect 
o\ = <7|| 7^ o~ j_ . Clearly, the noise terms could also have 
been expressed via £(x, t) = — V - ^(x, t), with correlations 
(C(x, i)C(x',i')> = -2A (cr||<9 2 + <J ± X7l) «5(x - x')tf(i - 0- 
Higher derivatives or powers of the order parameter which 
have been neglected above are expected to be small above 
the critical point. Near T c , they will be irrelevant, in the 
renormalization group sense, as we will see in Section 
4.2.1. 

This Langevin equation forms the basis of our theo- 
retical analysis. Since the isotropic diffusion term in Eqn 



(2.9), rV if, has been modified to (rud + T±X7j_)ip, with 



two different "diffusion" coefficients for the parallel and 
transverse subspaces, a serious question arises. Should we 
expect both, or just one, of these to vanish as T approaches 
criticality? For the equilibrium system, we recall that the 
lowering of r is a consequence of the presence of interpar- 
ticle interactions which slow (or even reverse) the decay of 
density gradients, as T decreases. In the two-temperature 
model, the parallel subspace is held at a higher tempera- 
ture than the transverse one, suggesting that density gra- 
dients should decay much faster in the parallel directions. 
Thus, we anticipate that, generically, rii > t±_ so that crit- 
icality, in particular, is marked by t± vanishing at positive 
711 . This is borne out by the structure of typical ordered 
configurations, namely, single strips aligned with the par- 
allel direction, indicating that "antidiffusion", i.e., t± < 0, 
dominates in the transverse directions below T r . 
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The inequality satisfied by the diffusion coefficients is 
obviously crucial for the correct description of criticality. 
Otherwise, it is fortunate t hat t he det ailed dependence 
of the parameters in Eqns ( [2.11 ) and ( 2.12 ) on the mi- 
croscopic J, T, an d T' , is n ot im portant. In the disor- 
dered phase, Eqns ( [2.11 ) and ( |2.12| ) can be simplified even 
further, by truncating the Landau expansion after linear 
terms, considering, in effect, a purely Gaussian theory. 

To con clude this section, let us check whether Eqns 
(2.11.2.12) satisfy the FDT. Comparing our noise corre- 
lations with Eqn (2.7), we read off N = ~o- a 8 a pd a dp. 
It is straightforward to see th at it is impossible to re- 
cast the right han d side of Eqn ( 2.11 ) in the Hamiltonian 
form of Eqn (2.8). We conclude that our continuum the- 
ory for the two-temperature model generically violates the 
FDT. However, in Section 4 we will show that the critical 
properties are controlled by a fixed point theory which is 
associated with, surprisingly, a system (not Ising) in equi- 



librium 19 . The implications are that FDT is restored in 
the critical region, for the leading singular parts of ther- 
modynamic quantities. 



3 Simulation Results 

This section will be devoted to the main results from 
Monte Carlo simulations on the two dimension model spec- 
ified in Section 2.1. In the first sub-section, we study the 
two-point correlations in the disordered phase, showing 
convincing evidence for power law decays in configuration 
space. In the second part, extensive investigations of the 
critical properties will be presented. We will use subscripts 
x and y for the "parallel" (||) and "transverse" (_L) direc- 
tions, respectively. In the remainder of this article, T will 
be quoted in units of the Onsager value T Q . 

3.1 Two-point Correlations Far Above Criticality 

For an Ising model in equilibrium, regardless of the dy- 
namical aspects of the system, two-point correlations 

G(X) EE ( SlSl+x ) 

display power law decay only at the critical point. How- 
ever, for a system driven into a steady state under non- 
equilibrium conditions, the asymptotic behavior of G will 
depend on the dynamics. In particular, for our model at 
any T above the critical point, G(x) decays as r -2 , where 
v ee |xj. Of course, the associated amplitude must be di- 
rection dependent, so that the integral over space is not 
divergent. In general, a simple rescaling of one co-ordinate 
(say, x) will bring it to the same form as an electrostatic 
potential of a quadrupole (in d = 2) J|,|]], so that there 
are long range positive/negative correlations in the longi- 
tudinal/transverse direction, i.e., 



G(x, 0) -> A/x 2 and G(0, y) -A'/y 2 



(3.1) 




Fig. 3.1. The structure factor 5(k) for a 128 x 128 lattice at 
T = 2.0. The integers (m x ,m y ) are wavenumbers: 128k/27r. 



difference (T' — T) as well, in a manner that vanishes as 
T" — » T. The presence of these amplitudes can be traced 
to FDT violation in these systems, as we will see in Sec- 
tion 4.1. Here, we provide some details of how this power 
law behavior can be observed directly in G, as well as the 
simpler, "indirect" method using the Fourier transform of 
G. 

We begin with the simpler way, i.e., by measuring the 
structure factor 



5(k)^]Te* kx G(x) 



and seeking a discontinuity singularity fll7[ at the origin 
k = 0. Unlike the Ornstein-Zernike form for the equilib- 
rium Ising model, where 

S eq (k)^ x + 0(k 2 ) 

with x being the static susceptibility, the limit depends 
on the direction of k as it approaches the origin ||[l7]] . In 
particular, the maximum range is embodied in, say, 



lim S(k x 



0, k y ) - lim S(k x ,k y 



0). 



(3.2) 



Another convenient way to express such a discontinuity is 



R EE 



lim |fc_ L Ho-5'(fc|| = 0,k±) 
limfcn_>o 6"(fe||, kj_ = 0) 



^1 



and A and A' are positive amplitudes depending on J and 
T. Had we used a finite T", both would depend on the 



Now, for a finite system, the wave-vectors k are discrete, 
being (k x , k y ) = 2n(jn x / L x , m y /L y ) with integer m*. Fur- 
ther, our system is set at M = 0, so that 5(0, 0) ee 0. Thus, 
through measuring 5(k) at the two lowest wave-vectors: 

(2n/L X7 0) and (0,2n/L y ) 

in a series of systems with various L x and L y , the discon- 
tinuity singularity in ( |3.2| ) can be estimated by extrapola- 
tion. Here, we are content to display the presence of such 
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Fig. 3.2. The longitudinal (o) and the transversal (•) corre- 
lation functions for a 256 x 256 lattice, at T = 3. To display 
data clearly, only x,y G [1,45] are shown. 



Fig. 3.4. Log-Log plot of 
a; -2 is also shown. 



-G(0, y) — constant. The asymptote 
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Fig. 3.3. Log-Log plot of G(x,0) — constant.The asymptote 
x~ 2 is also shown. 



a discontinuity, as in the first Monte Carlo study of driven 
systems [jl|. 

In Fig. 3.1, we show S(k) for the 128 x 128 case at 
T = 2.0 . For this plot, we made measurements every 
20 MCS on 5-10 runs, each up to 10 6 MCS. Note the 
sizable difference between S(2ir/L x , 0) and S(0,2ir/ L y ). 
By contrast, for the equilibrium case, these two quantities 
would be the same, within statistical errors. 

Parenthetically, let us point out that, for the m x = 
series, a multitude of long runs are crucial for minimizing 
the scatter in the data. This "difficulty" may be traced to 
anot her manifestation of long range correlations like Eqn 
(3.1), namely, the presence of strips along x which are 
anti-correlated in y. This property favors slow decays of 
strips with a range of widths, so that short runs may find 
only configurations "locked in" particular regions of phase 
space. We believe that this phenomenon is intimately tied 



to the metastability of "multistrip configurations" first re- 
ported in ]10| . 

Though the discontinuity in S(k) is the easiest way to 
deduce the presence of a power law decay, it is naturally 
more satisfying to observe the evidence directly. Indeed, 
in the early studies ^,|), G(x) is measured along the two 
axes and InG vs. lnr is plotted. However, the data used 
there to support the r~ 2 behavior are susceptible to criti- 
cisms. Typically, a straight line is drawn through less than 
10 scattered points, even though sizes up to L x = 300 
were used. In the study on the two-temperature model 
D, low statistics (20K MCS) further limits the reliability. 
To remedy this situation, we carried out a more careful 
study, measuring every 20 MCS on 5-10 runs, each up to 
10 6 MCS, on a 256 x 256 system at T = 3.0. In Fig. 3.2, 
we present G(x, 0) and G(0, y) vs. x or y. 

Apart from finding a more convincing evidence of the 
r~ 2 decay, we discovered a non-trivial finite size effect, 
which must be taken into account before a log-log plot 
can be exploited. Due to the conservation law M = 0, we 
have X) x G(x) oc M = 0, so that the asymptotic behavior 
of the two-point correlation in a finite system is not zero. 
Indeed, for the totally disordered state (T = oo), we must 
have G(x) = (N5 X , - 1) / (N - 1). As a result, for high 
T, the power law tail (3.1) can be easily overwhelmed by 
negative 0(1/ N) contributions. In practice, the precise 
values of these corrections (for finite T) are not known 
and represent the only parameters in the fitting procedure. 
The result reveals excellent r~ 2 decays. In Figs. 3.3 and 
3.4, dashed lines of slope —2 are included simply for the 
sake of compa riso n. From these plots, we conclude that 
the behaviors (JO) are well established. 



3.2 Critical Properties 

As T is lowered toward a critical value, T c , correlations 
build up and the system is expected to undergo a contin- 
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uous transition to an ordered state, characterized by phase 
segregation. Unlike in the equilibrium Ising case, this seg- 
regation occurs only along the y-direction, a phenomenon 
undoubtedly related to the anisotropic long-range correla- 
tion above T c . Another manifestation of this type of order- 
ing is that only 5(0, 2tt / ' L y ) diverges (as O(N) for T < T c ), 
while S(2ir/L x , 0) remains O(l) for all temperatures. The 
consequences of this behavior are quite serious. Further, 
without some theoretical elucidations, the pitfalls await- 
ing finite size analysis in computer simulation studies are 
well hidden. As an example, in contrast to the equilibrium 
case, the use of L x L samples (with a series of L's) will 
not lead to consistent collapse of data (Section 3.4). In- 
stead, due to the conserved dynamics and FDT violation, 
lengths / momenta in the two directions scale with different 
exponents ]l8| ], as will be shown in detail below (Section 
4.1). Thus, we resort to the method of "strong anisotropic 
scaling" i.e., using a series of rectangular samples 

with fixed L x /Ly + and non-trivial A. 

Before delving into the details of the results, let us 
provide an "appetizer" for the origins of a non-zero A, 
from "mean field theory". To describe this transition in 



the equilibrium case, Landau let r (in Eqn 2.4) g o th rough 



zero. As a result, the leading relaxation term in (2.4) van- 
ishes , re gardless of the noise "matrix" N. In particular, 
for ( |2.9[ ) both diffusion coefficients vanish at the same 
(critical) temperature. However, without the constraint of 
FDT, it is not necessary that both vanish at one T. Indeed, 
to describe phase segregation into horizontal strips (along 
with S(2tt/L x ,0) being 0(1)), it is nat ural to let only the 
transverse diffusion coefficient, t±, in (2.11), vanish. As a 



result, the leading restoration terms become —"S/^ip and 
d 2 tp, which translate into 

-dyip and d%(p, 

in the language relevant for our simulations. Since the ef- 
fects of these are supposedly comparable in the critical 
region, the implication is 

k x "~ ky . (3-3) 

But, the lowest wavevectors available to our finite system 
are (2ir/L x ,0) and (0,2tt / L y ). To insure that both types 
of fluctuations are comparable in the critical region, we 
must choose L x 
result 



Ly. Thus, we arrive at the "mean field" 



A 



MF 



1 



(3.4) 



As will be shown (Section 4.2), corrections to this result 
are expected to be small, even in d — 2 (unlike the uni- 
formly driven lattice gas, where A = 2). Thus, for con- 
venience, we simply used systems with A ~ 1, namely, 
20 x 20, 45 x 30, 80 x 40, 125 x 50, and 180 x 60. To have 
some confidence that 1 — A is indeed small and positive, 
we also carried out runs with 120 x 50 and 170 x 60. All 
results are entirely consistent within statistical errors, in 
that a good quality of data collapse is observed. By con- 
trast, the use of square samples does not lead to consistent 



data collapse (Section 3.2.3), signaling the presence 
other scaling variable: 



L x /L 



l+A 

y 



of an- 



(3.5) 



Finally, let us remark on possible choices for order pa- 
rameters. Due to the conservation law, M is fixed at zero 
and cannot serve as the order parameter. Instead, in the 
seminal work on driven systems [Q and several subsequent 
studies a quantity was used which is sensitive to or- 
dering into horizontal strips rather than vertical ones: 



E 



E 



(3.6) 



However, (3J3) is likely to overestimate the segregation 
into two macroscopic phases |pd| , |l3| . For example, it can- 
not distinguish single strip configurations from multiple 
strip ones, as long as the latter are all horizontal. Further, 
since it is a sum of S(k) over all k's (along both axes in k 
space) up to the ultraviolet cut-off, its scaling properties 
are quite nebulous. Therefore, we decided upon the struc- 
ture factor S(0,2Tr/L y ), being the clearest measure of the 
spontaneous symmetry breaking, as our order parameter 
Jul . Configurations with multiple strips contribute little 
to this quantity. Also, as a correlation function associated 
with the smallest possible wavevector, there is little doubt 
on how to compare simulation data with field theoretic 
predictions. 



3.2.1 Crossing of a Cumulant Ratio, Data Collapse and the 
Exponent v 

For the equilibrium Ising model, Binder [^H introduced 

a convenient ratio g(L,T) = 3 — (^(X) s ) 4 ) / ^E s ) 2 ) i 
which varies monotonically between 2 and 0, as T ranges 
in [0, oo], regardless of L. In the thermodynamic limit, 
it should be a simple step function, with the step occur- 
ring at the critical temperature. For finite L's, g(L,T) 
curves should cross at temperatures approaching T c . Fur- 
thermore, the value at the crossing is predicted to be uni- 
versal pi[ |, so that this method has been the favorite for 
the first determination of the critical temperature. For our 
model, s i s clearly futile, since it is identically zero. In- 
stead, we must consider the first non-zero wavenumber in 
the Fourier transform: 



E< 



with k = (0, 2n/Ly) . 



Note that we chose ^|s| J / (L x L y ) as our order parame- 
ter. The equivalent cumulant ratio here is pd|,^3[ 



9(L V ,T) 



(3.7) 



In this definition, the dependence on L x is suppressed, 
since we have used only sizes which correspond to fixed 
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a ee L x /L y +A . Also, due to s* (k) = s(-k) ^ s(k), it 
is more convenient to have 2 as the first term, so that 
g(L y , oo) = 0. At the opposite end, g(L y ,0) — 1, so 
that g(L y — > oo, T) is precisely the Heaviside function: 
0(T c — T). In the inset of Fig. 3.5, we plot the ratio against 
T for the system sizes chosen. Within the errors of our 
simulations, the crossing occurs at T c ~ 1.37. 
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Fig. 3.5. Scaling plot of g(L y ,T) vs. t=L 1/v (T - T c )/T c , for 
different system sizes L x xL y : | 180x60, x 170x60, * 125x50, 
□ 120x50, ■ 80x40, o 45x30, • 20x20. Inset: g vs. T. 



More importantly, when plotted against the scaling 
variable 



/ EE I ^ ~ 1 ) 



(3.8) 



all curves may be collapsed into a single one (Fig. 3.5), 
with 



and 



T c ~ 1.370(2) . 



i/~ 0.62(3). 



(3.9) 



(3.10) 



Note that we have tested the effects of having A slightly 
less than 1 as well. We are encouraged by the data asso- 
ciated with these cases (120 x 50 and 170 x 60) also lying 
within the typical scatter of the "universal" curve. Thus, it 
is not crucial that we keep careful account of the possible 
presence of a small (1 — A). The value for the exponent v 
is entirely consistent with the results of renormalization 
group calculations (Section 4.2.3). 

Finally, we emphasize that a general scaling ansatz for 
g should be of the form g(L x , L V ,T) — g(a,t). Thus, the 
use of square samples would have introduced the added 
complication of a varying a. However, if a were to enter 
g as an additional power, then data collapse can still be 
achieved, but with another value for v. We will return to 
this topic in Section 3.2.3. 



3.2.2 Scaling Plot of 5(0, 2tt/L v ) and the Exponent r\ 

From scaling arguments, the other independent exponent 
may be chosen as r\. By definition, it is the anomalous 
dimension associated with the structure factor at criti- 
cality in the limit of small k. In models with "strong 
anisotropy" such as ours, we must be more careful and 
define it to be the anomalous dimension in 5(0, k y — > 
0;T = T c ) -► ky 2+r i. (See Section 4.2.2 below for other 
77-like exponents.) For that reason, 5(0, 2n/L y ; T) serves 
conveniently in a finite size scaling analysis to extract 77. 
In practice, it is more convenient to use O] 



S = S(0,2ir/L y ;T)^- S m 2 (j- 



which is normalized to unity at T = 0. Thus, the presence 
of scaling is signaled by SL x L y L~ 2+v ~ SL y +n being 
a function of i alone. Of course, there are generally two 
branches, associated with T's above/below criticality. For 
T > T c , we may define 7 in the usual way, via S(T) — > 
(T — T c )~ 7 in the thermodynamic limit. In other words, 
SL A+V should be proportional to i~ 7 . Using the definition 
of t, we obtain the usual scaling relation 2 — r; = j/u. 
Meanwhile, for T < T c , 5 plays the important role of 
^M 2 ) for non-conserved systems and serves to probe /3, 

the order parameter exponent. So, SL^f^ u should be a 
function of t alone, leading to another scaling relation: 
A + r] = 2(3 jv. Thus, data collapse in a single scaling 
plot, for both T above and below T c , is a good test of the 
hyperscaling relation: 



2/3 + 7 = v(2 + A) . 



(3.11) 



As will be shown below, field theory predicts A — l — rj/2, 

so that its validity can be tested by plotting SLl +v ^ 2 
against t and checking for data collapse. Further consis- 
tency can be verified by extracting the exponents 7 and (3 
from the two branches and checking if they correspond to 
(2 — r\)v and (1 + r]/2)v/2 1 respectively. 

Fig. 3.6 shows a log-lo g plo t of the bes t dat a collapse 
of the two branches, using (3.9) for T c and ( 3.10 ) for v. To 
guide the eye, and to show consistency, we have inserted 
lines with the expected slopes. From this graph (and oth- 
ers with varying values of rf) we estimate the value 



77 ~ 0.13(4) 



(3.12) 



and conclude that a consistent set of exponents exists, 
including 

/3~ 0.33(2), and 7- 1.16(6). 

The smaller samples deviate from the collapse as the tem- 
perature rises. This deviation can probably be traced to 
the following. For small system sizes, larger temperature 
ranges are needed to accumulate data for the scaling plot, 
corresponding to T's further outside the critical region. 
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Fig. 3.6. Scaling plot of SL y +v ^ 2 against t for different system 
sizes L x xL y : | 180x60, x 170x60, * 125x50, □ 120x50, ■ 
80x40, o 45x30. 



Fig. 3.7. Scaling plot of g(LxL,T) vs. (T/T c - l)L 1/l/ , with 
T c =1.34 and ^=1.00 for different system sizes L x xL y : □ 
80x80, I 60x60, x 40x40, * 20x20. 



3.2.3 Data from Square Samples 

In order to test the importance of "strong anisotropic 
scaling," we have also performed simulations using square 
samples from 20 x 20 to 80 x 80. With data from such 
runs, it is possible to achieve reasonably good collapse of 
g(L x L; T), using T c ~ 1.34 and v ~ 1.00 (Fig. 3.7). How- 
ever, with these values for T c and v, the best "collapse" we 
can achieve for S is shown in Fig. 3.8, with r\ ~ 0. Better 
collapse of these data can be accomplished, but only by 
using values of T c and v different from those in Fig. 3.7, or 
at the expense of a large, negative value of rf. We believe 
that the lack of good data collapse (with a consistent set 
of values for the critical temperature and exponents) is a 
strong indication that another scaling variable is playing 
a dominant role. In this case, it would be L y /L 1 x +A . If we 
use the mean- field value for Zi, then this variable ranges 
over a factor of 4, for the set of data presented. Such a 
large variation should be enough to explain the significant 
scatter in Fig. 3.8. 

In the literature, it has been claimed that consistent 
data collapse can be achieved using square samples ||. 
However, we b elie ve that, by using a different order pa- 
rameter (Eqn (3.6)), these authors may be probing a dif- 
ferent "phase." In connection with the uniformly driven 
lattice gas, the possible existence of such a "new" phase 
has been conjectured recently |32| and may resolve the dif- 
ferences observed here. Further investigations are clearly 
necessary to clarify these issues. 



3.2.4 Fluctuations and Specific Heat 
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Fig. 3.8. Scaling plot of SL y +rl ^ 2 against i for different system 
sizes L x xL y : □ 80x80, | 60x60, x 40x40, * 20x20. 



different temperatures. Nevertheless, we may continue to 
regard (H) as an internal energy, associated with the in- 
terparticle interactions. Furthermore, since anisotropy is 
expected, we may define the "longitudinal" and "trans- 
verse" parts of this energy, by {Tt x ) and (H y ) respectively, 
where 



T~La 



-47 



n(i x ,i y )n(i x + l,i y ) 



and 



-47 n (i x , i y ) n (i x , i v + 1) 



(3.13) 



For an Ising model in thermal equilibrium, both the free 
energy and the internal energy, (H) , are well defined and 
useful concepts. For non-equilibrium systems in steady 
state, the concept of free energy becomes less clear, es- 
pecially in our case of having thermal reservoirs at two 



Needless to say, these are nothing but the energies stored 
in the broken bonds in the two directions. 

As in equilibrium, the internal energy is a function of 
T and we may define its response to changes in T as the 
heat capacity C 
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Fig. 3.9. Fluctuations in TL y , the "transverse" part of the Fig. 3.10. Log-log plots of the heights of the peaks (from Fig. 
energy, vs. 
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T, for different system sizes: | 180x60, x 170x60, 
• 120x50, □ 80x40, ■ 45x30, A 20x20. 



3.9) vs. L y . Inset shows ln(l — T p /T c ) vs lnL y , with T p being 
the positions of the peaks. The line is of slope —1/0.62. 



C 



d(H) 
dT 



Similarly, its fluctuations can be measured. To emphasize 
the similarities and differences between our system and 
one in equilibrium, we consider 



C F = x 



((Any 



k B T 2 ■ 

Replacing TL in these expressions by H. x and Ti y , we can 
define similar quantities: C Xl C Fx , etc. Of course, had we 
used T" = T, we would have C = C F = 2C X = 2C Fx , etc. 
Thus, measurements of all these quantities will provide not 
only a specific-heat-like exponent a, but also the degree 
of FDT-violation. 

We are able to compile reasonably good statistics on 
the fluctuations (C^'s). For the critical region, the best 
data are associated with C Fy , the fluctuations in the "trans- 
verse" part of the energy (Fig. 3.9). By contrast, due to 
the positive long range correlations, there are fewer broken 
"longitudinal" bonds, which may be the reason behind the 
lower quality of the data for C Fx . To study scaling prop- 
erties, we consider first the positions of the peaks, T p (L y ). 
In Fig. 3.10 (inset), we show a log-log plot of 1 — T p /T c vs. 
L y , along with a line of slope of —1/0.62, so that we have 
a consistent estimate for v. On the other hand, a simi- 
lar plot for the heights of the peaks (Fig. 3.10) leads us 
to a/u ~ 0.28. This value of a is also entirely consistent 
with the scaling relation predicted from field theory (Sec- 
tion 4.2.3), i.e., 2 — a = (d + A) v. Unfortunately, scaling 
plots for C Fy do not exhibit data collapse of the same qual- 
ity as those for g and S. We believe that better statistics 
would remedy this situation. Similarly, the data for energy 
are not precise enough for estimates of the heat capacities 
(C's). Nevertheless, we may conclude from these observa- 
tions that the qualitative features of the C's are the same 



as the Cf's. In particular, a log- log plot of the heights of 
the C' y maxima vs. L y provides us with a similar a. How- 
ever, it is not surprising that the numerical values of the 
C's are quite distinct from those for the fluctuations, since 
FDT is clearly violated. 

We end this section with a brief comment on energy 
fluxes. Since our system is coupled to two temperature 
baths, we expect that there would be a steady flow of en- 
ergy through our system, from the bath with T" to the one 
with T. To be specific, with each spin-exchange update, 
we can track whether it is associated with an x-bond or 
a y-bond. Respectively, the energy change resulting from 
these exchanges would be determined by T 1 and T. We 
have verified that, in the steady state, the average energy 
change associated with an x-bond spin-exchange is posi- 
tive. Indeed, far from being just fluctuations, the average 
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Fig. 3.11. Energy flux {8Ti, x } — (87i y ), vs. T, for various system 
sizes: + 180 x 60, x 125 x 50, * 80 x 40, □ 45 x 30. 
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change in a MCS, which we denote by (6H X ), is a few per- 
cent of 1(H) , for all the sizes we investigated. Of course, 
being in steady state, the average change associated with 
y-bond exchanges is of the same magnitude, but nega- 
tive: (SH y ) = — (SH X ). In other words, the energy flux 
through the system, (SH X ) — (SH y ), is an extensive quan- 
tity. By contrast, for an equilibrium system, the flux is 
fluctuating around zero, at the order of y/N. Similar dif- 
ferences have been observed in uniformly driven diffusive 
systems [[33]. Qualitatively, this flux rises with tempera- 
ture until about the critical point. Thereafter, it appears 
relatively constant (Fig. 3.11). A quantitative study of this 
flux, which lies outside the scope of this paper, would un- 
doubtedly yield rich information about an important non- 
equilibrium characteristic of this system. 



4 Field Theoretic Studies 

Exact solutions for many-particle systems far away from 
thermal equilibrium are restricted mostly to one-dimen- 
sional cases Q with excluded volume interactions alone. 
Thus, progress for higher-dimensional models relies en- 
tirely on mesoscopic continuum theories, such as the Lange- 
vin equation for our two-temperature model. While these 
field theories cannot be rigorously derived from the mi- 
croscopic dynamics, many of their key ingredients and 
predictions are easily tested by Monte Carlo simulations. 
Considerable confidence in these theories can therefore be 
established before they are used to predict less easily mea- 
surable quantities. In the following, we will outline how 
some of the simulations results of the preceding sections 
are reflected in our continuum theory. We begin with the 
disordered phase, before turning to fully developed critical 
behavior. 



4.1 Long-Range Correlations 

Well above T c , typical configurations are nearly homoge- 
neous, and fluctuations of the local magnetization </?(x, t) 
away from zero are very small. It is therefore reasonable to 
neglect fourth-order gradient terms as well as higher pow- 
ers of if in our Langevin equation. The resulting theory, 
for the equilibrium Ising model, is completely analytic for 
all T > T c . In contrast, for our far-from-equilibrium sys- 
tem, even the disordered phase already exhibits nontrivial 
singularities. In the following, we will discuss how these 
manifest themselves in our continuum theory. 



We begin with the simplified version of Eqns ( 2.11 2.12 ) 
which is fully appropriate for the disordered phase: 

dt<p(x, t) = \{t ± VI^ + T {l d 2 ip} - V ■ ry(x, f) (4.1) 

The noise correlations are given by 



(7 ?Q (x,i)^(x',i')> = 2\<r a S a p6{x-x')5(t-t') 



(4.2) 



Since Eqn (4.1) is linear, it is easily solved via a Fourier 
transform to momentum and frequency space. With <p(x, t) 



J p(k, t) exp*(k ■ x + ut), where J .= / A I % 



we find 



-ik • ?7(k, lu) 



ico + xfrxkl+T^ +0(fc4,fc4 jfc 2 fc 2) 



(4.3) 

The corrections in the denominator are a reminder of the 
neglected fourth-order derivatives. Arbitrary correlation 
functions follow from (ETq) by averaging appropriate prod- 
ucts of <p(k, oj) over theuaussian distribution of the noise. 

The most remarkable property of the disordered phase 
is the presence of power-law correlations, even well above 
T c . In the equal-time structure factor, being the Fourier 
transform of the correlation function, these should be re- 
flected as an anomaly near the origin. Starting from the 
full dynamic structure factor, defined via (y(k, Lo)ip(k' , a/)) 
= S(oj + Lu')S(k + k')iSo(k, u>), its Fourier transform back 
into the time domain, So(k, t), is easily found: 



So(M) 



M-MMk,o)) 



a±k 2 A 



T ± kl+T l{ k 2 + 0(kl,kl,klkl) 



x exp 



-A(r±fci+T||fcf) \t\ 



(4.4) 



The subscript reminds us that we are considering a purely 
Gaussian theory. Setting t — yields the equal-time struc- 
ture factor, S'o(k) = So(k, 0): 



So(k) 



en kl +ankf. 



Tiki 



0{k i j fetj ^ k | II 



(4.5) 



In contrast to what one might have expected naively, (4.E ) 
is not a simple anisotropic generalization, such as S'e(k) oc 
l/(r + k\ + bk 2 ), of the usual isotropic Ornstein-Zernike 

structure factor l/(r+fc 2 ) for the equilibrium Ising model. 
The key differ enc e resides in the discontinuity singularity 
exhibited by (4.5). Defining 



R = 



_ lim|fc_L|- 



+0 S(fc|| =0,k x ) 



limk||_,o<S , (fe||ikL = 0) 



we find that both S^k) and its isotropic version lead to 
R = 1, indicating that such an anisotropy is too wea k to 
affect the long- wavelength limit. In stark contrast, (4.E) 
results in R = (o~ut±)/(o~±t\\), since the value of Sq at the 
origin depends on the angle under which k vanishes! 

At a more fundamental level, this discontinuity singu- 
larity is a direct consequence of the violation of the FDT. 
For linear theories such as ours, the FDT demands that 
the matrices of noise correlations, o~ a 8 a f3, and of diffusion 
coefficients, T a 5 a p, should be proportional to one another. 
Thus, if the FDT were to hold, the general form (4.5) 
would be immediately reduced to 5 e (k). This is the case 
for, e.g., an Ising model with anisotropic interactions. For 
our two-temperature model, however, there is no such con- 
straint, so that generically R ^ 1. At criticality, the FDT 
is maximally violated, i.e., R diverges as tj_ — > 0. 
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One consequence of such a singularity in Sq is that 
its Fourier transform, the two-point correlation function 
G(r), becomes long-ranged, decaying as r~ 2 at large dis- 
tances. The amplitude, in addition to being proportional 
to (R — 1), has a dipolar angular dependence, resulting 
in positive G(xu,x± = 0) and negative G(xu = 0, Xj_), in 
agreement with our data. Of course, our continuum the- 
ory presupposes an infinite system, so that no finite size 
corrections appear here. Moreover, the dipolar character 
combined with the 0(fc 4 ) corrections ensures that an ap- 
propriate angular average of G is again exponentially cut 
off However, the divergence of R at criticality indi- 
cates a crossover to a different behavior, which will be the 
subject of the next section. 

To conclude, we note that the Ising "up-down" symme- 
try of our model ensures that all three-point functions van- 
ish identically above and at criticality. This feature clearly 
distinguishes our model from the uniformly driven system 
where such correlations are generically nonzero and can 
even be singular at the origin, along selected directions 



4.2 Critical Properties and the Nontrivial Fixed Point 

As the temperature approaches T c , the local magnetiza- 
tion begins to develop large fluctuations, both in abso- 
lute value and in gradie nts. This signals the breakdown of 
the linear theory, Eqns ( |4.l| , [4.2[ ), as a valid description of 
the disordered phase. We sho uld t herefo re ret urn to the 
full Langevin equation, Eqns ( 2.11 ) and ( 2.1 2| ) . Our first 
task is to identify all couplings relevant for critical prop- 
erties, i.e., retaining non-zero values at the fixed point of 
the renormalization group. We will see that simple dimen- 
sional analysis alone can already exhibit some of the key 
differences between our two-temperature system and its 
equilibrium counterpart, Model B. Most remarkably, the 
fixed point theory restores the FDT, but with respect to 
a different Hamiltonian. We define a set of critical expo- 
nents and perform a scaling analysis for our model, by 
anticipating the scaling forms for a few key quantities. Fi- 
nally, we provide the technical foundation for this scaling 
analysis, in the form of an explicit renormalization group 
calculation, up to and including two-loop order |19|. 



4.2.1 The Fixed Point Theory and FDT Restoration. 



In the continuum theory, (2.11,2.12), criticality is marked 
by the vanishing of the transverse diffusion coefficient, r±, 
while its parallel counterpart tn remains positive. Recon- 
sidering our Landau expansion in powers of gradients, we 
clearly need the (V^) 2 ^ term to limit fluctuations with 
large transverse gradients, since T±'S/ 2 L ip cannot serve this 
purpose near T c . In contrast, fluctuations with large par- 
allel gradients are still controlled by T^d 2 ip, so that the 
higher order term, (d 2 ) 2 ip, is not needed and will be ne- 
glected. As a result, the two leading linear terms in the 
Langevin equation, near criticality, are (V^) 2 ^ and T\\d 2 tp. 



Introducing a characteristic momentum scale for trans- 
verse critical fluctuations, |kj_| ~ /i, we conclude, on the 
basis of simple dimensional analysis, that parallel wave 
vectors scale as |fcii | ~ |k^| 2 ~ /j, 2 . In the long-wavelength 
limit (/i — > 0) which is the focus of our study, parallel 
gradients are therefore less relevant than transverse ones: 
d 2 ip 3 may be neglected in favor of V 2 ^ 3 , and transverse 
noise correlations dominate parallel ones. Similarly, the 
mixed term q x 9 2 V^ is irrelevant compared to (V^) 2 ^. 
Collecting so far, we arrive at a reduced Langevin equa- 
tion, containing relevant terms only: 



d t <p(x,t) = A{(r ± - V 2 L )V 2 L ^^ 
+C(x,i) 



3! 



(4.6) 



Here, the noise has been restricted to the transverse sub- 
space, with correlations 



(C(x, <)C(x', f )) = -2AtT ± V 2 L ^(x - x')6(t - t') 



(4.7) 



This equation forms the starting point for the analysis of 
universal critical properties. Continuing the dimensional 
analysis, we find that Xt ~ /i~ 4 , characteristic for a con- 
served order parameter. Care must be exercised with the 
scaling of the spatial volume. For example, S(x — x') ~ 
since parallel lengths generate an extra factor /it -1 . 
By construction, both rn and a± are of order 1. In fact, 
the latter is a trivial coefficient which will be rescaled to 
1 in the following. The role of tii is more intriguing and 
will be discussed in Section 4.2.3. It is now easy to ob- 
tain £(x, t) ~ ^( d + 7 )/ 2 so that the order parameter itself 
scales as ip ~ /i^ -1 )/ 2 . The upper critical dimension of 
our model is d c — 3 since u ~ /i 3 ~ d . Since all irrelevant 
couplings have been neglected, we refer to Eqns (4.6,4~7j), 
somewhat loosely, as the fixed point theory. Strictly speak- 
ing, we should set t± = also, since the naive scaling 
ti ~ [i 2 indicates that zero is the fixed point value of the 
transverse diffusion coefficient. 

It is instructive to define D = d + n, where n is the di- 
mension of the parallel subspace. In the spirit of the more 
general models alluded to in Section 2.2, n may be con- 
sidered a model parameter. Of course, we require d > n, 
for a transverse subspace to exist. Then, the magnetiza- 
tion scales as p,( D ~ 2 ^ 2 and u ~ /i 4_jD , whence D c = 4. 
In this form, the similarity to Model B scaling, where 
n = 0, is quite apparent. However, we emphasize that the 
physical upper critical dimension is d c — D c — n, below 
the Ising model value 4. This already indicates that the 
two-temperature model falls outside the Ising universality 
class. Clearly, the only "interesting" case is n = 1: then, 
nontrivial exponents are expected in two dimensions, as 
borne out by our Monte Carlo simulations. In contrast, 
the choice n = 2 results in an upper critical dimension of 
2, for a model that is defined only in d > 3 dimensions. As 
a result, only mean-field exponents should be observable 
in this case. 

Returning to Eqns ( [l.6|j4.7|) , we observe an intriguing 
consequence of neglecting irrelevant terms: the reduced 
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Langevin equation for the near-critical theory obeys the 
FDT. Defining N = — AV^, and noting that N is a positive 
definite symmetric operator, we can rewrite our Langevin 
equation as follows: 



%>(M) = -N^ + C(x,t) 



(4.8) 



with a Hamiltonian that is most easily expressed in mo- 
mentum space: 



HM = 



/> ( - k) 



k\ 



^(ki)^(k 2 )^(k 3 )^(k4) 



ii 
"4! 



'ki,k2,k 3 ,k4 

x(2TT) d 6(k 1 + ... + k 4 ) 



(4.9) 



Since the only relevant contribution to the noise acts in 
the transverse subspace, its correlations are simply 



(C(x,t)C(x',0)=2M(x-x')5(t-t') 



(4.10) 



This structure leads to several important consequences. 
First of all, near criticality, the static properties of the 
system are in fact equilibrium-like, being controlled by 
the distribution exp(— Ti). Interestingly, even though our 
microscopic dynamics is purely local, effective long-range 
interactions are generated at the mesoscopic level. These 
are dipolar in nature, as reflected by the static Gaussian 
propagator 



5 (k) = 



ki 



HI 



(4.11) 



In fact, it is helpful to recall the Hamiltonian for uni- 
axial ferromagnets with dipolar interactions J36|. Assum- 
ing that the spins are aligned with the parallel direction 
and located at r and r', the dipolar interaction takes 



dx 2 . 



/r 



d+2 



the form s r U(r — r')s r ', with U(r) oc 

A Fourier transform to momentum space yields t/(k) cx 
k 2 /k 2 . Including the exchange interaction, we obtain the 

propagator for uniaxial dipolar ferromagnets, Sd(k) oc 
1 -l 

+ k 2 + bk 2 jk 2 , where r — > marks the critical point 

and b is a constant. At first sight, the k-dependences of 
these propagators, S<j(k) and So(k), are obviously not 
identical. However, closer examination reveals that the dif- 
ference resides in irrelevant terms alone (e.g., fcjj, in a sum 

with kj_)\ Of course, the nonlinear term, in both cases, is 
just the usual ip A interaction. Thus, we conclude that the 
leading static critical singularities of our two-temperature 
model fall into the universality class of an equilibrium 
Hamiltonian, which also describes uniaxial dipolar ferro- 
magnets. This result is quite remarkable, given that we 
started from an FDT- violating microscopic dynamics. At 
a technical level, it simplifies our renormalization group 
analysis, since the statics of uniaxial dipolar ferromagnets 
has been discussed in the literature 



Second, the FDT generates a hierarchy of equations, 
relating correlation and response functions. These are more 
easily discussed if we first recast our theory in terms of 
a dynamic functional. The latter also represents a much 
more convenient starting point for the RG calculation 
of the dynamic critical properties. Introducing a Martin- 
Siggia-Rosc |3£| response field y>(x, i), we obtain the dy- 
namic functional at the fixed point, 



J[(f>M 



dtd d x 



6T~L 1 

(pdtf + <pN- pNip } 

Sip J 



(4.12) 



This form has the advantage that both correlation and re- 
sponse functions can be computed as functional averages 
with weight exp(— J). Specifically, response functions are 
just expectation values containing the field <p(x,t). For 
example, (ip(x' ,t')ip(x,t)) represents the response of the 
local magnetization tp, at (x, t), to an infinitesimal per- 
turbation at (x',t'). If we were to add a magnetic field 
to the Hamiltonian 7i, the dynamic zero-field susceptibil- 
ity would be just x(k, to) = k\ (<p(— k, — u>)p(k, u/)) where 
the prefactor reminds us that the order parameter is con- 
served and therefore cannot respond to uniform perturba- 
tions (such as a homogeneous magnetic field). By virtue 
of the FDT, this susceptibility is related to the dynamic 
structure factor, according to 



S(k,w) = 2XUJ- 1 im[x(k,w)] 



(4.13) 



Higher order response and correlation functions satisfy 
similar relations. 

Some ca ution must be exercised here. The FDT, and 
hence Eqn (4.13), holds strictly only at the fixed point of 
the two-temperature model. Away from the fixed point, 
the full functional is generically FDT- violating, including 
irrelevant operator s whi ch cannot be recast in the form 
prescribed by Eqn (4.12). Thus, these irrelevant operators 
generate different corrections-to-scaling for response and 
correlation functi ons. W ith respect to simulation data, 
this implies that (4.13) holds only inside the critical re- 
gion, for the leading singular parts of S and \- A similar 
statement qualifies the relation between energy fluctua- 
tions and the specific heat. In this sense, "FDT breaking" 
operators are dangerously irrelevant: their absence mimics 
a symmetry which is not strictly satisfied. 

Before concluding this section, let us briefly return to 
the discussion of more general models and consider the 
case of more than two temperatures, T\ > T2 > ... > T<£. 
Naturally, we expect the lowest of these to control criti- 
cality. Let us assume that, say, the m lowest temperatures 
are degenerate and are decreased while keeping the re- 
maining n = d — m (higher) temperatures fixed. In the 
continuum limit, this will give rise to a series of diffusion 
coefficients, t\ > r 2 > ...t„ > r n+ i = ... = Trf. Not surpris- 
ingly, r„ + i, ...,Td will vanish first, defining the transverse 
subspace. The remaining n — d—m dimensions specify the 
parallel subspace, marked by finite diffusion coefficients at 
criticality. 
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4.2.2 Scaling Analysis 

Near criticality, large fluctuations on all length scales dom- 
inate the behavior of the system, so that renormaliza- 
tion group techniques are indispensable. These allow us 
to compute, e.g., scaling forms and critical exponents ex- 
plicitly, in an expansion about the upper critical dimen- 
sion. Here, we will anticipate the scaling properties of our 
model, deferring technical details to the next section. 
The dis cus sion leading to the fixed point theory, Eqns 



scale as |kj_| ~ 
obey fcii - l 1+A 



tV' , whence parallel momenta 



(4.6) and (4.7), suggests that the critical behavior of the 
two-temperature model is distinct from its Ising origins: 
since parallel and transverse wave vectors scale with differ- 
ent powers, the upper critical dimension is shifted to d c = 
3. Anticipating renormalization, we reformulate this scal- 
ing as | foil | ~ |kjj 1+zl , introducing the strong anisotropy 
exponent A. For comparison, weak anisotropy refers to 
systems where A remains zero and only scaling ampli- 
tudes are anisotropic. Equilibrium examples for the for- 
mer include, e.g., Lifshitz points or structural phase tran- 
sitions H^J , while the Ising model with anisotropic interac- 
tions falls into the second category. Far from equilibrium, 
the usual driven lattice gas Mis the prototype model for 
strongly anisotropic scaling |ll| , |l3| . Since all characteris- 
tic lengths are expected to scale with the same exponents 
near criticality, a finite size scaling analysis should rely on 
aspect ratios with constant L\i/L_^~ A . 

For any system with strong anisotropy, irrespective of 
its universality class, the renormalization group predicts 
the general scaling form of, e.g., the dynamic structure 
factor near criticality: 

S(M;tx) = l~ 2+ " n S{kJl 1+A ,'kx/l,tl z ;Tjl 1 / v ) (4.14) 



Here, I is just a scaling factor. Eqn ( 4.14 ) can be viewed 
as a definition of the critical exponents v, z, r\ and A. The 
latter is a new exponent, in addition to the usual two 
independent static exponents v and 77, and the dynamic 
exponent z. Different universality classes are distinguished 
by the characteristic values of these exponents, expressed, 
e.g., through their e-expansions. These will be discussed 
in the next section. 

The scaling of the structure factor is particularly in- 
structive, since strong anisotropy plays a key role here. 
For example, four 77-like exponents can be defined @, 
all of which would take identical values in an isotropic 
system. In contrast, here we should introduce r\\_ and 

via S(fc|| = 0,k_L,£ = 0;rj_ = 0) ~ |k x |~ 2+,u , and 
S(k {l ,k ± = Q,t = 0;r± = 0) ~ |fc|| |- 2+ "ii , for |k| -► 0. 
Similarly, the two-point correlation function decays with 
exponents n' ± and 77^ for large distances, G(r\\ , rj_ = 0, t = 

0;tl = 0) ~ r ~ < - d ~ 2+v 'w\ an d G(r\\ = Q,rj_,t = 0;t± = 



0) - r {{ y T ,lJ . Fortunately, even though these four ex- 
ponents are generically different, they are not indepen- 
dent. Instead, they are related to A and 77 through simple 
scaling laws j|. For example, rjx =77, and n'« — 71 fj^ 3 - ) . 
Similarly, we can define two critical exponents, v±_ and 
mi: choosing I to satisfy t±/1 1 / u — 1, transverse momenta 



Tj_' ' = r_|_ , resulting in i/± = v and 
z/|| =i/(l + A). The scaling of the momenta with time 
t determines two dynamic critical exponents, z±_ and z« , 
according to |kj_| ~ t~ x / z± and k\\ ~ t _1 /*u. We easily 
find z±_ = z and z\\ = z/(l + A). The dynamic exponents 
are especially interesting here, since they have not been 
investigated previously, unlike the static ones which agree 
with those for the dipolar system |3?J . 

To compare with simulation data, several other expo- 
nents are interesting. The order parameter exponent (3, 
obtained from an equation of state, controls the response 
of the magnetization. In the absence of dangerous irrele- 
vant operators, it can be extracted from a scaling analysis, 
discussed in the next section: 

P=\v{d-l + ^) 

Other exponents of interest include the "specific heat" 
exponent a. To be specific, let us consider fluctuations of 
the energy, Cf, first. Taking the naive continuum limit of 
Ji, the average total internal energy would be identified 
as (/ d d x </3 2 (x, t)) (plus higher order, less singular terms). 
Thus, the fluctuations are associated with the connected 
part of J d d xd d y (<p 2 (x, t) <p 2 (y, t)). To extract the singu- 
lar part of this operator, we only need to consider the fixed 
point functional. But the latter corresponds to a Hamil- 
tonian system, so that the computation follows standard 
routes for static systems. In this sense, there is no need 
to recompute a, since it will be related to v through the 
usual scaling relation. The only point we should emphasize 
is that, due to the anomalous scaling of the longitudinal 
momenta ~ |kjj 1+zi ~ /i 1+A ), we have 



2 - a = (d + A) v 



(4.15) 



Since FDT is restored for the fixed point, we may also 
conclude that the singular part of the heat capacity, C, 
is of the same form, so that there will not be a new ex- 
ponent. Of course, in contrast to an equilibrium system, 
C is not strictly equal to Cf, so that different amplitudes 
associated with the leading singularities may be expected. 

We may further inquire into the singular parts of the 
energy "difference" (H x ) — (H y ) in simulations, which is 
identically zero in case the two temperatures are equal. 
In higher dimensions, the equivalent quantity would be 
(7i±) — (d — 1) \H\\). Taking the naive continuum limit, 
the corresponding operator is 

(/ {(V^) 2 -(d-l)(cV) 2 }^ 

With a higher naive dimension than the total energy, its 
singular parts should be of higher order, so that verifying 
their presence would not be facile.. From the theoretical 
perspective, perhaps a more interesting question is to ex- 
plore the energy flux: (STi x ) — (5H y ), which is associated 
inherently with a non-equilibrium process. However, since 
the fixed point dynamic functional is Hamiltonian, we be- 
lieve that the singular parts of such operators would be 
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(a) 



(b) 



(O 



(d) 



Fig. 4.1. Diagrammatic representation of the bare propagator 
(a), correlator (b), vertex (c), and insertion (d). 



irrelevant (though dangerously so). Though intriguing, the 
tasks of deriving and analyzing their corresponding oper- 
ators in field theory lies outside the scope of the present 
paper. 



4.2.3 Renormalization Group Computation at Two-loop 
Order 

In this Section, we finally provide some technical details of 
the renormalization gro up anal ysis. The st arting point is 
the Langevin equation (16 4/7), recast as ( 4.12 ), the dy- 
namic functional p6| . Following standard methods E(J , 
we perform a renormalized perturbation expansion, orga- 
nized in powers of the nonlinearity. We first collect the 
elements of perturbation theory, i.e., the bare correlation 
and response propagators and the vertex. We then focus 
on the one-particle irreducible vertex functions and iden- 
tify those which are primitively divergent. Using dimen- 
sional regularization, we compute the associated poles in 
e = d c — d, up to and including two loops, followed by 
the renormalization of vertex functions and coupling con- 
stants. Finally, we establish and solve the renormalization 
group equation for the vertex functions, thus arriving at 
the full scaling behavior. 

The bare propagators are easily read off from the Gaus- 
sian theory, corresponding to u = 0, indicated by the sub- 
script (-) . For the response propagator, shown in Fig. 
4.1a, we find 



G (k,w) ee (<p(k,u)cp(-k,-(j)) 

1 



ILU + X ( Tx_k\ + T||fcjj + kj_\ 



(4.16) 



while the correlation propagator (Fig. 4.1b), also referred 
to as the correlator, is given by: 



C (k,w) = (tp{k,uj)tp(-k, -w)) 
2Xk 2 , 



iuj + X (r±_k\ + T\\k 2 + kj^ 



(4.17) 



The latter, of course, is just the dynamic struct ure f actor 
for the Gaussian theory. A comparison of Eqns ( 4.14 ) and 



(4.17) allows us to identify the Gaussian, or mean-field, 
values of the critical exponents: 



v = — , z = 4, 
2 



i] = and A = 1 



(4.18) 



Clearly, these will acquire nontrivial corrections of 0(e). 

The theory possesses only a single nonlinearity, eas- 
ily identified from Eqn ( 4.12; ) as Xu J dtd d x(V 2 L (p)(p 3 . In 
Fourier space, this gives rise to a four-point vertex, shown 
in Fig. 4.1c, 



—Xuk i 



(4.19) 



in Feynman integrals. Here, k± is the (transverse) mo- 
mentum carried by the <^-leg. However, the expansion is 
not organized simply in powers of u, as one can infer from 
an additional scaling symmetry of the theory: Due to the 
spatial anisotropy, parallel and transverse momenta (or 
lengths) can be rescaled independently, namely fen — > fci i /a 
while kx remains unchanged. The functional, Eqn ( 4.12 ), 
is invariant under this transformation, provided we rescale 



'II 



a 2 Tii and u — > au, so that ut. 



-1/2 



is the invariant 



form of the nonlinearity. For convenience, we also absorb 
a geometric factor 4ir into the definition of the effective 
coupling constant: 



1 

IT UT \\ 

47T II 



-1/2 



(4.20) 



In the following, we consider the critical theory, i.e., we 
set tj_ to zero in all propagators. However, t± does require 
renormalization which can be extracted from insertions 
of Xt± J dtd d x('V 2 L ifi)(p into graphs of the critical theory. 
As usual, the insertion corresponds to a two-point vertex, 
shown in Fig. 4. Id, carrying nonzero external momentum 
For later reference, we note that vertex functions with 
insertions can be resummed, to generate the vertex func- 
tions of the theory above T c Q . 

Even though the detailed momentum dependence of 
propagators, vertex and insertion is of course different 
from Model B, the topology and combinatorics of the 
Feynman diagrams are identical. This helps to simplify 
some of the technicalities. For example, denoting bare one- 
particle irreducible vertex functions with N (N) external 
tp- (tp-) legs and L insertions as N . L , it is straightfor- 
ward to identify the primitively divergent ones as A,i;0> 
A.3:0> and A,i ; i- At the upper critical dimension, i.e., in 
d = 3, the last two of these both scale as fc^. This mo- 
mentum dependence is already carried by a factor of fc^ 
on one of the external legs, so that the remaining integrals 
are momentum-independent. In contrast, A,i;0 has the di- 
mension of k\ , so that the integral itself must contribute 
another factor of k\ , in addition to the one carried by the 
external leg. Similar to Model B, each of the three vertex 
functions gives rise to one nontrivial renormalization, so 
that we anticipate two independent exponents provided an 
infrared stable fixed point can be found. All other vertex 
functions, including especially /2,0;0; are primitively con- 
vergent. It is now straightforward to evaluate these dia- 
grams and extract the singularities. Relegating the details 
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to the Appendix, we only quote the results 

1 



A,i ; o(k±,0) 

A.a.o^ki.Q}) 
A,i;i(k±,0;k /5 0) 



Afcl 



1 



u 2 kT 2e J 



-\uk 1 



(4.21 



3 

--uh 
-\T ± k 2 ± \ 1 - -ul x 



-u 2 I 2 



-u 2 I 2 



3u 2 I 2 



u 2 h 



Here, J as well as I\ and I2 are integrals defined in the 
Appendix. 

Next, we renormalize our theory, by introducing a set 
of renormalized quantities (which are identified by the su- 
perscript (R)) and the associated Z- factors. These are re- 
lated to the original (bare) set of quantities by: 



T ll = z v lz n T \ 



(R) 



= Z-'Z^r 



(R) 



A = (ZjZs^ZxXW 



(4.22) 



/j. £ Z Z u u 



(R) 



Note that is dimensionless, and can be used to define 
a renormalized effective coupling 



,.(«) 



-in 



XR) 



(4.23) 



which is, unlike u, also dimensionless. Based on (p( R ' and 
ip( R \ we consider a set of renormalized vertex functions: 



11 

= Z^ 2 Z§ /2 r^ N . L ({k,u;};X,T h T^u). (4.24) 

The renormalization conditions consist of demanding that 
all Ij -ty's b e fin ite as e — * 0. As a result, through Eqns 
( 4.21 ) and ( 4.22| ), all renormalized quantities are implic- 
itly dependent on /1 - the (momentum) scale at which 
these conditions are imposed. Of the many routes to real- 
ize finiteness, the minimal subtraction scheme is the most 
facile. Only pure poles in e need to be introduced in the 
Z-factors, to cancel those in the -T's. Thus, we can easily 
read off the results for the Z's. 

First, we find to all orders in e that 



Z^Z^p 



1 



(4.25) 



i.e., the couplings A and rii require no renormalization 
while the Z-factors for the fields are inverse to one an- 
other. The remaining Z-factors must be determined order 
by order in perturbation theory. To two loops, we obtain 



Z^ = 1 



Z T± = 1 



Z u = 1 



-J9 2 
6 y 

1 - 

nhg- 



\hg 



0(g 3 ) 

4 Jl _ 2 h 
15 p f 

T 1 ~ 32 



0(g 3 ) (4.26) 
-0(g 3 ) 



where the expressions I\ and I2 are defined in the Ap- 
pendix. 

\ Once the (finite) I^^'s are known, we can find their 
scaling behavior in the infrared limit by studying their 
"flow" as fi — > 0. Since the bare -T's are independent of /i, 
flow equations Q for the f^'s arise from applying fid^, 
at fixed bare quantities, to Eqn (4.24): 



= fid 



V\ bare r « { + (3(g)d g + ...} (4.27) 

where (3(g) = ^d^\ bare g and the dots represent similar 
terms for the other parameters. Defering the complete 
equation until below, we focus on (3(g), which controls 
the flow of the effective coupling g. To describe critical 
behavior, we seek infrared stable fixed points of the flow, 
i.e., zeroes of (3(g) which are attractive under fi — > 0. In 
our case, 

m = '9 - \g + \Kg 2 + 0(.g 3 )| (4.28) 

where K = 17/54 + In (2/V3). For d > 3, the only stable 
fixed point is the Gaussian g* = 0. It becomes unstable 
below the critical dimension, i.e., in d = 3 — e, to be re- 
placed by a second, nontrivial fixed point which emerges 
smoothly from the origin: 



(4.29) 



(4.30) 



(4.31) 



l + -Ke + 0(e 2 ) 

Near g* , g(fi) flows as 

ffO) = g* + [s(mo) - g*] O/mo)" 

if fi = fiQ initially. Here, 

u = e- -Ke 2 + 0(e 3 ) 
3 

is positive (within perturbation theory), so that ( 4.29| ) is 
indeed stable. This is also a critical exponent, controlling 
corrections-to-scaling . 

Having established the existence of an infrared stable 
fixed point below the upper critical dimension, we now 
focus on the scaling properties of our theory which can be 
extracted from the RG equation. The latter expresses the 
//-independence of the bare vertex functions, at fixed bare 
parameters, as a differential equation for the renormalized 
vertex function. To begin with, let us introduce several 
additional Wilson functions: 



7. 



fid^Z, = (3(g)d g In Z m and k, = fid^ ln( 



(4.32) 



where the symbol • represents any of the variables if, (p, 
T||, ti, or A. Given the definitions of the couplings and 
the Z-factors computed above, these functions satisfy the 
following identities, valid to all orders in perturbation the- 
ory: 



7A = 7r,| = , 7 V = -7^ , K T „ = -«A = l<p , 



(4.33) 
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Thus, there are only two independent Wilson functions, and 
7 V and k Tj _, whose expansions are easily computed up to 
0(g 2 ). For the following, only their values at the fixed 
point are important: 



7 *^( 5 *) = - e 2 + 0( £ 3 ) 

3 

54 " 



< *N 1 2 

*tA9 ) = n e + n 



K 



o(e 3 ; 

-,(R) 



(4.34) 



of the 



Finally, we recall that the vertex functions _T:~ , 

-" N,N-L 

critical theory, with insertions of t± , are easily resummed 
to give the vertex functions rifl T of the disordered phase, 

at finite tj_ > 0. Suppressing the superscript (i?) on all 
parameters for the sake of clarity, we quote the full RG 
equation: 



= 



j/xc^ + f3{g)d g + txk Tj _8 t± + T||/« T|| 9 T|| + \n\d\ 

N N } 

~ y7 v j^v,Jv({ k ' w }; r |h T -i-- A,.g,/i) (4.35) 



This equation can be solved, using the method of char- 
acteristics. In the immediate vicinity of the fixed point, 
where the theory exhibits scaling, we can incorporate di- 
mensional analysis and the anisotropic scale transforma- 
tion into the solution, whence: 



r jv,Ar({ k M; 



(4.36) 



= Fr^ N ({Vi ± l-\k ll l- 2+ ^/ 2 ,ujl- i+ ^};Txl- 2+K ') 



Here, the overall scaling exponent is 



p = d + 5 



N 



(d 



1 + ^7*) 



N 



{d + 3 



(4.37) 



and t he param eter I is just a scaling factor, as in Eqn 
( 4.14 ). In Eqn ( 4.36 ), we list only those parameters in the 
argument that take an active part in the scaling form. 

Let us illustrate how Eqn (4.14) emerges from the so- 
lution of the RG equation. To obtain the scaling form of 
the structure factor, we recall the connection between the 
two-point correlations and the two-point vertex functions: 



S(k,uj;T ± ) 



r 2 , (k,a;,Tx) 
|A,i(k,w; t±) 



Inserting the appropriate scaling forms, and performing a 
Fourier transform from the time- into the frequency do- 
main yields Eqn (4.14): 



S(M;t x ) = 

= r 2 +^*s(k ± r \ fc|, r 2 +^2^4- 7 * . Tir 2+ K * } 

In analogy with the scaling form of Model B, we now define 
two static exponents 



V = 1 



(4.38) 



1 


2 


— K* 


1 


1 


2 


+ 12 



1 

36 



67 

54 



In- 



+ 0(e 3 ) (4.39) 



which are the only two independent indices. We emphasize 
that both of these differ from their counterparts for the 
Ising model near d = 4. Thus, the two-temperature model 
falls into a different universality class, defined by the dipo- 
lar system as shown by a comparison with the exponents 
quoted in Hj. However, a two-loop calculation is neces- 
sary in order to exhibit the differences between Model B 
and the dipolar system. At the order of one loop, rj retains 
its mean-field value, and v is determined by combinatorics 
alone. Moving beyond statics, we find that the dynamic 
exponent z is related to the static ones by the conservation 
law: 



' 7 



While this scaling law is also obeyed by the Model B ex- 
ponents, the explicit result for z is characteristic for the 
two-temperature model, as well as the dipolar system in 
the presence of a conservation law for the magnetization. 
Finally, we quote the strong anisotropy exponent A, which 
controls the anisotropic scaling of parallel and transverse 
momenta: 



A 



i 1 * 



1 

1-2*7 



Scaling laws relate all other exponents to these two, e.g., 
the "specific heat" exponent a as discussed above. 

The most interesting remaining index is the order pa- 
rameter exponent /3 which can be extracted from the equa- 
tion of state. Since there are no dangerous irrelevant op- 
erators here, a simple scaling analysis is sufficient. In or- 
der to describe the effect of a magnetic field, we add a 
term A J dtd d x<fW 2 L h(x±) to the dynamic functional. The 
transverse gradient operator imposes the (relevant) effect 
of the conservation law. Clearly, this will generate a non- 
vanishing, spatially varying magnetization, M(k), which 
breaks the Ising symmetry of our system. Introducing 
the vertex generating functional r{tp,<p}, the equation 
of state for the renormalized quantities follows from 



d S 
dk\ bCp 



-^{^^}|k=0,^=0, V =M 



Here, the conservation law, reflected in the derivative djdk\, 
selects the first nonzero Fourier component of the right 
hand side. Expanding r in powers of M, the right hand 
side is reduced to a sum of vertex functions which are 
computed in the disordered phase: 



dk 2 



M N 

N 



fen =0,o; = 0};ti 



The scaling form of the equation of state now follows di- 
rectly from the solution of the RG equation, Eqn ( 4.36 ): 



h(M,T^) = i( d + 3 -hy^ h {Mi-( d - 1+ ^y 2 ,Txi- 1/v ) 
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Choosing I such that Ml ^ d 1 +l r i)/ 2 — 1 ) the equation of 
state is recast in standard form 

h(M,T X ) = M 5 f(T X M- 1/!3 ) 

with an appropriately defined scaling function /. We can 
now read off the magnetic-field exponent 

rf + 3-f 

and the order parameter exponent 

(3 = \v{d-\ + \). 



5 Summary and Outlook 

We have presented a detailed study of a non-equilibrium 
Ising system, with Kawasaki exchange dynamics coupled 
anisotropically to two thermal baths. Even when one of 
the baths is set at infinite temperature, this system is 
known to display a second order phase transition similar to 
the equilibrium model, with generic scale invariance at all 
temperatures above criticality p|,^2[. Both field theoretic 
renormalization group techniques and Monte Carlo simu- 
lations on 2-dimensional square lattices have been used. 
Excellent agreement between the two approaches is found 
for T > T c . Since this system displays strongly anisotropic 
scaling properties, we exploit finite size scaling with rect- 
angular samples to control the effects of a ne w s caling 
variable associated with the aspect ratio (Eqn 3.5). The 
results are entirely consistent with the predictions of the 
field theory. By contrast, data from square samples fail 
to collapse, showing that the aspect ratio indeed plays a 
significant role. 

One of the more intriguing aspects of this system is 
the following. The fixed point controlling critical proper- 
ties can be associated with a uniaxial magnet with dipolar 
interactions in equilibrium. More specifically, the leading 
singularities in the thermodynamic functions of our non- 
equilibrium lattice gas are, in the critical region, identical 
to those for the equilibrium magnetic system. Of course, 
our lattice gas is still a bona-fide non-equilibrium system. 
To find these aspects, we must study either sub-leading 
singularities, such as corrections to scaling, or quantities 
associated with dangerous irrelevant operators [142] . A good 
example of the latter is the energy flux through the sys- 
tem, which is easily measured in simulations but quite 
difficult to compute theoretically. 

Looking beyond this study, we see opportunities for 
further pursuit. Apart from the obvious need for better 
statistics and larger system sizes, there are many new av- 
enues. To end this paper, we list a few. (a) Very little 
is known about this system below the critical point. Only 
the fluctuations of the interface were scrutinized |Q . Cor- 
relations within each bulk phase are expected to be long- 
ranged, but confirmations from simulations are still lack- 
ing. Similar to the case of a uniformly driven lattice gas, 



domain growth during phase segregation is expected to 
display strongly anisotropic time scales H|. Unlike that 
case, there is no particle-hole asymmetry in our model. 
As a consequence, there may be no serious disagreement 
between the microscopic lattice gas model and the meso- 
scopic Cahn-Hilliard like approach, (b) For the driven lat- 
tice gas, imposing open j45| or skewed-periodic (4(| bound- 
ary conditions leads to surprisingly different behavior, es- 
pecially at low temperatures. Viable field theories for these 
systems are yet to be formulated, (c) On the other hand, a 
complete renormalization group analysis exists for multi- 
temperature models, as shown above. Though most of the 
predictions are mean-field like, it would be interesting to 
carry out simulations to confirm or disprove them. 

Of course, this model represents only a small exam- 
ple in the vast field of systems in non-equilibrium steady 
states. Other models subjected to two thermal baths abound 
E7j , not to mention complex physical systems ranging 
from heating a pot on a stove to the ecosystem of the 
earth. Hopefully, the richness displayed in this very simple 
two-temperature model will spur further interest in pur- 
suing other models and, most importantly, enhance our 
insights into this vast class of physical systems in non- 
equilibrium steady states. 
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6 Appendix 



In this appendix, we include some details of how to extract 
the singular parts of the three vertex functions of inter- 
est: A,i;0j A,3;0j and Figs. 6.2-4 show the corre- 
sponding Feynman graphs, up to and including two loops. 
Contributions that vanish in dimensional regularization 
or due to causality have already been omitted. We stress 
that we are only interested in the ultraviolet singular parts 
of the vertex functions, in the critical theory. To regular- 
ize the infrared behavior of our integrals, it is sufficient 
to keep transverse external momenta finite. Parallel ex- 
ternal momenta and external frequencies may eventually 
be set to zero, since they contribute to finite terms only. 
The abbreviated notation f = f f ^ will be used 
throughout. 

The structure of A,i;0 at this order is: 
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A,i ; o(k, Si) = iS2 + \ (ki + T {l kf) -^{Xu) 2 k 2 ± x 
qlj_ Go (qi , w a ) C (q 2 , uj 2 ) C (q 3 , uj 3 ) 



X(27r) d+1 5(^ qi -k)^< 



J?) 



(6.1) 



Any neglected terms involve at least three loops. This ex- 
pression can be simplified considerably. First, it is straight- 
forward to integrate over oj's, and the external S2 may 
be set to zero. Next, we absorb the parallel diffusion co- 
efficient Tii into the parallel mom enta, so that the effec- 
tive coupling constant, Eqn ( 4.20 ), appears explicitly. Fi- 
nally, we can recast the remaining momentum integrals in 
a more symmetric form, resulting in: 



A,i ; o(k,0) = A Ik 



— Au fcjJ(k) 
6 



(6.2) 



where the remaining integral (which still carries a nontriv- 
ial k-dependence, as we will see below) is defined as 



J(k) 



qi,--,q3 



<?3_L 



3 



(?2X + «2||J ('/',_ + <R 

5 (qi)5o(q 2 )5o(q 3 )(2 7 r) d ( 5(^q l -k).(6.3) 
q1.q2.q3 i = i 

In the second line, we have recognized that the momenta 
under the integral combine into t hree factors of the static 
Gaussian structure factor, Eqn (4.11). Thus, the result- 
ing expression, up to a prefactor AfcJ , is identical to the 
static two-point vertex function of the uniaxial dipolar fer- 
romagnet, at two-loop order, with the correct momentum 
integrals and combinatoric factors! This structure reflects 
the restoration of the FDT, Eqn ( 4.1 3| ) , and can be ob- 
served at each order in perturbation theory 

Similar reductions occur for the other two vertex func- 
tions. Let us demonstrate it for the one-loop graph of the 
(dynamic) four-point function, -Ti,3;0- Its calculation is fur- 
ther simplified by choosing transverse external momenta 
at a symmetry point, defined by 

k l± ■ k j± = k 2 h± (4% - 1) /3 , 
for i,j — 1, ...4. Up to one loop, the expression reads: 
r i, 3i o({ki, S2i = 0}) = -\ukl x + 3(Xu) 2 kl ± x 
x/ qf± Go(qi,wi)5o(q2,w 2 ) x 

J qi,...,ui2 

x (27r) d+1 5(qi + q 2 - ki - k 2 )<5(wi + lo 2 ) 
Defining 

q.2 ki 



(6.4) 



/i({k ( ,0}) 



q 2 ± q 2 ± (27r) d £( qi 



<?i||) (<?2X + <? 2 ||) 
5o(qi)5 (q2)(27r) d 5( qi + q 2 - k x - k 2 ) (6.5) 



Fig. 6.1. Feynman diagram contributing to A,i;0- 



we write 



A,3; ({k l ,0}) 



-Xuk 1J _ + —Xuitk 1± Ii (6-6) 



Again, we recognize the zero-frequency part as the one- 
loop form of the static four-point function (up to a factor 
Xk 2 ± ), involving just static propagators. The explicit re- 
ductions at the two-loop level involve combinations of sev- 
eral dynamic diagrams into one static one. For brevity, we 
only quote the results, so that the leading singular parts 
of Pi 3 read: 



A, 3; o({k 4 , 0}) = -Xukl ± \ 1 - \uh + lu 2 h 2 + M 2 I 2 



where the last integral is given by 



/ 2 ({k i) 0}) 



qi,-,q4 i= i 



(6.7) 



(6.8) 



x (27r)^<5(qi + q 2 - ki - k 2 )c5(q 3 + q 4 - qi - k 3 ) 

Finally, we turn to the last vertex function, namely 
A,i;i which carries the insertion. As in Model B, the in- 
tegrals contributing to A,i ; i are just those of A,3 ; o : The 
insertion corresponds to two of the four external legs be- 
ing "tied together" . In general, it carries both momentum 
(kf) and frequency (wr), but the latter can be set to zero 
again. Thus, 

ri,i.i(k,0;kj,0) = -Xr ± k 2 ± jl - X -uh + \u 2 ll + hi 2 I 2 

(6.9) 

It is now straightforward, if tedious, to evaluate these 
integrals and to extract the singularities. The one-loop 
diagrams present no difficulties. It is rather instructive, 
however, to consider one of the nontrivial two-loop dia- 
grams as an example. We illustrate our procedure with 
the help of the integral contributing to /~i,i;0 : 

3 

J(k) = / (2K) d 5(J2 q* - k ) So(qi)So(q 2 )So(q 3 ) . 

"'q1.q2.q3 i—i 

(6.10) 



'qi,q2 



First, we integrate over the parallel momenta. To do so, 
we recast the 5-function for the parallel momenta via its 
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XX XXX 



xb X) 




x>- xX>-- 

Q a 



Fig. 6.3. Feynman diagrams contributing to 



Fig. 6.2. Feynman diagrams contributing to A,3 ; o- 



Fourier transform, leading to a product of three integrals 
of the form 



/ 



— 7 q± V exp(^g||) = iexp(-giH) (6.11) 



2tt 



The integral over x is now easily performed, and the 
external ku may again be set to zero, resulting in 

* 1 f 1 
J (kj_,0) = - / 7-5 5 5— r X 

3 

x(2^) d -M(^q 4 x-k ± ) (6.12) 

i=l 

The advantage is obvious: the extremely inconvenient q\_- 
terms in the denominator have been reduced to quadratic 
(gj_) powers, in the process cancelling the (^-factors in 
the numerator. The remaining integral is of a much sim- 
pler form and requires only the familiar tools for evaluat- 
ing standard </> 4 -theory integrals |4C| . In a similar fashion, 
one can simplify the two-loop integrals for fi,3 ; o- Collect- 
ing our results for the integrals, we obtain, up to O(e) 
corrections: 

1 



J(k ± ,0) = - — (ki) 1 -' 



II ({]*,<)}) = 1 {l + \ [-C E + M127T)]} (klj 



i-e/2 



/2({ki,0}) 



i-C B +m(8V37r) | (A* 



Here, Ce = 0.577... is Euler's constant. The anticipated 
factor of kj_ is, finally, displayed explicitly in J(kj_, 0). It 
indicates that this inte gral is a correction to the zero-loop 
term \k\ in Eqn ( |6.2| ): 



A,i ; o(kx,0) = Afci--A{ t 2 A:2 
6 



36e 



(fc 



2 \l-e 



[l + 0(e) 



1 _ Lu 2 k- 2t j 
6 - 1 



In the minimal subtraction scheme, only the pure e- 
poles in the I's and J will be needed. Thus, we define 



J 
h 



if 



1 

~36e 
1 

27 

^{l + £ [-C £ 

i l n 

8? 



1 



+ ln(127r)]} 
+ ln(8\/37r) 



(6.14) 



(6.13) 



where all O(l) corrections have been neglected. Note that, 
within this framework, I\ is "more than the square of" 1\ . 
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